
#define _USE_MATH_DEFINES
#include "math.h"
#include "malloc.h"

class airpotentials
{


public:

	// Math functions
	double scalprod(double a[3], double b[3])
	{
		double res = 0;
		for(int i=0; i<3; i++) res+= a[i]*b[i];
		return res;
	}

	void vectprod(double a[3], double b[3], double *res)
	{
		res[0] = a[1]*b[2]-a[2]*b[1];
		res[1] = a[2]*b[0]-a[0]*b[2];
		res[2] = a[0]*b[1]-a[1]*b[0];
	}


	airpotentials()
	{

	}

	// Model parameters
	double Alpha, Beta, R0;
	int ssize[3][4];
	int terms[4];
	int terms_total;
	int terms_max;
	double *K;

	double *set_term[4];
	double *set_term_D[4];

	void LoadModelFromFile(char *modelfile)
	{
		
		FILE *mf;
		mf=fopen(modelfile,"r");

		fscanf(mf,"%lf %lf %lf",&Alpha,&Beta,&R0);

		for(int j=0;j<4;j++)
		{
			terms[j] = 0;
		}

		for(int i=0;i<3;i++){
			for(int j=0;j<4;j++)
			{
				int n;
				fscanf(mf,"%d",&n);
				ssize[i][j] = n;
				terms[j] += n;
			}
		}
		
		terms_total = terms[0]*terms[1]*terms[2]*terms[3] + 2;

		terms_max = 0;
		for(int j=0;j<4;j++)
		{
			if (terms_max < terms[j]) terms_max = terms[j];
		}

		K = (double*)malloc(terms_total*sizeof(double));

		for(int j=0; j<4; j++)
		{
			set_term[j] = (double*)malloc(terms_max*sizeof(double));
			set_term_D[j] = (double*)malloc(terms_max*sizeof(double));
		}

		for(int i=0;i<terms_total;i++){
			fscanf(mf,"%lf",&(K[i]));
		}
		fclose(mf);
	}
	

	double ComputeEnergy(double x[4][3], double *forces) // U in K, forces in N*10^-10
	{

		const double epsilon_const = 0.0001;

		//////////////////////////////////////////////////////////////////////////
		// 1. Angles and distances computation
		//////////////////////////////////////////////////////////////////////////

		double r[4][3];
		double rs[4];
		for (int a = 0; a < 3; a++)
		{
			r[0][a] = (x[2][a] - x[0][a]);
			r[1][a] = (x[3][a] - x[0][a]);
			r[2][a] = (x[2][a] - x[1][a]);
			r[3][a] = (x[3][a] - x[1][a]);
		}

		for(int i=0;i<4;i++) rs[i] = sqrt(scalprod(r[i],r[i]));

		// Molecules centers
		double xc[2][3];
		for (int a = 0; a < 3; a++)
		{
			xc[0][a] = (x[0][a] + x[1][a])/2.0;
			xc[1][a] = (x[2][a] + x[3][a])/2.0;
		}

        // Intermolecular distance
		double Rv[3]; // intermolecular vector
		double Re[3]; // normalized Rv

        for (int a = 0; a < 3; a++)
		{
			Rv[a] = (xc[1][a] - xc[0][a]);
		}
		
		double R = sqrt(scalprod(Rv,Rv));

        for (int a = 0; a < 3; a++)
		{
			Re[a] = Rv[a]/R;
		}

		// Angles
		double gamma[3];
		double ax[4][3]; // atoms vectors (relative)
        for (int a = 0; a < 3; a++)
		{
			ax[0][a] = x[0][a] - xc[0][a];
			ax[2][a] = x[2][a] - xc[1][a];
		}

		double l[2] = {sqrt(scalprod(ax[0],ax[0])),sqrt(scalprod(ax[2],ax[2]))};

		for(int j=0; j<3;j++){ 
			ax[0][j] = ax[0][j]/l[0];
			ax[1][j] = -ax[0][j];
			ax[2][j] = ax[2][j]/l[1];
			ax[3][j] = -ax[2][j];
		}

		int an[2] = {0,2}; // atom number for gammas
		double sc0 = scalprod(ax[0],Re);
		double sc1 = scalprod(ax[2],Re);
		if (sc0>=0) { an[0] = 0; } else {an[0] = 1; sc0 = -sc0;}
		if (sc1<=0) { an[1] = 2; sc1=-sc1; } else {an[1] = 3;}

		if (sc0>1) sc0=1;
		if (sc1>1) sc1=1;
		gamma[2] = acos(sc0);
		gamma[1] = acos(sc1);

		double n0[3],n1[3],tao1[3],tao2[3],y1[3],y2[3],eta1[3],eta2[3];

		// normals to planes
		vectprod(Re,ax[an[0]],n0);
		vectprod(Re,ax[an[1]],n1);

		double n0s = sqrt(scalprod(n0,n0));
		double n1s = sqrt(scalprod(n1,n1));

		int g0_defined = 1, n0_defined = 1, n1_defined = 1;

		if (n0s > epsilon_const) 
		{
			for(int j=0;j<3;j++) n0[j] = n0[j]/n0s;
		} else
		{
			g0_defined = 0;
			n0_defined = 0;
		}

		if (n1s > epsilon_const) 
		{
			for(int j=0;j<3;j++) n1[j] = n1[j]/n1s;
		} else
		{
			g0_defined = 0;
			n1_defined = 0;
		}

		if (g0_defined == 0) 
		{
			// gamma_0 isnt determined, let gamma_0 = 0
			gamma[0] = 0.0;
		} else { 
			double cosn1n2 = scalprod(n0,n1);
			if (cosn1n2>1) cosn1n2 = 1;
			if (cosn1n2<-1) cosn1n2 = -1;
			gamma[0] =  acos(cosn1n2) - M_PI/2.0;
		}

		//////////////////////////////////////////////////////////////////////////
		// 2. Energy computation
		//////////////////////////////////////////////////////////////////////////

		double expr[4], expr2[4];
		
		for(int i=0;i<4;i++) expr[i] = exp(-Alpha * rs[i]);

		// U_SR terms
		double S[2] = {0,0};

		for(int ri=0;ri<4;ri++)
		{
			expr2[ri] = expr[ri]*expr[ri];
			S[0] += expr[ri];
			S[1] += expr2[ri];
		}

		// U_LR terms
		for (int anglenum = 0; anglenum < 4; anglenum++)
		{
			int num = 0;
			for (int angle = 1; angle <= ssize[0][anglenum]; angle++)
			{
				set_term[anglenum][num] = 1;
				set_term_D[anglenum][num] = 0;
				num++;
			}

			for (int angle = 1; angle <= ssize[1][anglenum]; angle++)
			{
				if (anglenum < 3)
				{
					set_term[anglenum][num] = cos(angle * gamma[anglenum]);
					set_term_D[anglenum][num] = -angle *sin(angle * gamma[anglenum]);
				}
				else
				{
					set_term[anglenum][num] = pow(R, -angle);
					set_term_D[anglenum][num] = -angle *pow(R, -(angle+1));
				}
				num++;
			}

			for (int angle = 1; angle <= ssize[2][anglenum]; angle++)
			{
				if (anglenum < 3)
				{
					set_term[anglenum][num] = sin(angle * gamma[anglenum]);
					set_term_D[anglenum][num] = angle *cos(angle * gamma[anglenum]);
					num++;
				}
			}

		}

		int n = 2; 
		 // Potential energy (PE)
        double PE = 0, PEfrom2 = 0;
		double force_gammaR[4] = {0,0,0,0};

		for (int q1 = 0; q1 < terms[0]; q1++)
			for (int q2 = 0; q2 < terms[1]; q2++)
			{
				double pC01 =    set_term[0][q1] *   set_term[1][q2];
				double pC0D1 = set_term_D[0][q1] *   set_term[1][q2];
				double pC01D =   set_term[0][q1] * set_term_D[1][q2];

				for (int q3 = 0; q3 < terms[2]; q3++)
					for (int q4 = 0; q4 < terms[3]; q4++)
					{
						
						double pC23 = K[n] * set_term[2][q3] * set_term[3][q4];

						PEfrom2 +=  pC01 * pC23;
						force_gammaR[0] += pC0D1 * pC23;
						force_gammaR[1] += pC01D * pC23;
						force_gammaR[2] += K[n] * pC01 * set_term_D[2][q3] * set_term[3][q4];
						force_gammaR[3] += K[n] * pC01 * set_term[2][q3] * set_term_D[3][q4];

						n++;
					}

			}
     

		double Rfact = pow(R/R0,Beta)*exp(-Beta*R/R0); //pow(R,Beta);

		PEfrom2 *= Rfact;
		for(int ang=0;ang<4;ang++) force_gammaR[ang] *= Rfact;

		PE = PEfrom2 + K[0]*S[0]+K[1]*S[1];



		///////////////////////////////////////////////////////////////////
		// 4. Forces computation
		///////////////////////////////////////////////////////////////////

		vectprod(n0,ax[an[0]],tao1); // tao1/l0
		vectprod(n1,ax[an[1]],tao2); // tao2/l1

		vectprod(n0,Rv,y1); // y1
		vectprod(n1,Rv,y2); // y2

		vectprod(n0,ax[an[1]],eta1); // eta1/l1
		vectprod(n1,ax[an[0]],eta2); // eta2/l0

		double sing[3] = {sin(gamma[0]),sin(gamma[1]),sin(gamma[2])};
		double cosg[3] = {cos(gamma[0]),cos(gamma[1]),cos(gamma[2])};

		double CF1 = 0;
		double CF2 = 0;
		double CF3 = 0;
		double CF4 = 0;

        if (g0_defined == 1) {
			double n0m = n0s*l[0]*R;
			double n1m = n1s*l[1]*R;
			
			for (int j = 0; j < 3; j++)
			{
				tao1[j] *= sing[0]*l[0]/n0m;
				eta2[j] *= l[0]/n0m;
			}

			for (int j = 0; j < 3; j++)
			{
				tao2[j] *= sing[0]*l[1]/n1m;
				eta1[j] *= l[1]/n1m;
			}		
		
			CF1 = sing[0]/n0m;
			CF2 = sing[0]/n1m;
			CF3 = 1.0/n0m;
			CF4 = 1.0/n1m;
		}


		double CF5 = cosg[1]/R;
		double CF6 = cosg[1]/l[1];
		double CF7 = cosg[2]/R;
		double CF8 = cosg[2]/l[0];



		// U partial derivatives
		double force_r[4];
		double force_R;

		for(int ri=0;ri<4;ri++)	force_r[ri] = -Alpha*(K[0]*expr[ri] + 2*K[1]*expr2[ri]);

		force_R = Beta*(1.0/R - 1.0/R0)*PEfrom2 + force_gammaR[3];

		double Ufact[3] = {0,0,0};
		if (cosg[0]>epsilon_const) Ufact[0] = - force_gammaR[0]/cosg[0]/2.0;
		if (n1_defined == 1) Ufact[1] =   force_gammaR[1]/sing[1]/2.0;
		if (n0_defined == 1) Ufact[2] =   force_gammaR[2]/sing[2]/2.0;

		// Forces to atoms

		double force_atom[4][3];
        for (int j = 0; j < 3; j++)
        {
			double ffR =  Re[j]*force_R/2.0; // R force and Morse
			double ffr20 = r[0][j]*force_r[0]/rs[0]; // 2-0 force
			double ffr30 = r[1][j]*force_r[1]/rs[1]; // 3-0 force
			double ffr21 = r[2][j]*force_r[2]/rs[2]; // 2-1 force
			double ffr31 = r[3][j]*force_r[3]/rs[3]; // 3-1 force

			double ffgamma0[4] = {0,0,0,0}; // Gamma0 forces
			if (g0_defined == 1) {
				double LCF1 = eta1[j] + eta2[j];
				double LCF2 = tao1[j] + tao2[j];
				double LCF3 = CF1*y1[j]   + CF3*y2[j];
				double LCF4 = CF4*y1[j]   + CF2*y2[j];

				ffgamma0[0] = Ufact[0] * ( - LCF1 - LCF3 - LCF2); 
				ffgamma0[1] = Ufact[0] * ( - LCF1 + LCF3 - LCF2); 
				ffgamma0[2] = Ufact[0] * (   LCF1 - LCF4 + LCF2); 
				ffgamma0[3] = Ufact[0] * (   LCF1 + LCF4 + LCF2);

			}

			double ffgamma1[4] = {0,0,0,0}; // Gamma1 forces
			if (n1_defined == 1) {
				double ffgamma1_C1 = Ufact[1] * (  ax[an[1]][j]/R + CF5*Re[j] );
				double ffgamma1_C2 = Ufact[1] * (   CF6*ax[an[1]][j] + Re[j]/l[1] );

				ffgamma1[0] = ffgamma1_C1; 
				ffgamma1[1] = ffgamma1_C1; 
				ffgamma1[2] = - ffgamma1_C1 - ffgamma1_C2; 
				ffgamma1[3] = - ffgamma1_C1 + ffgamma1_C2; 
			}

			double ffgamma2[4] = {0,0,0,0}; // Gamma2 forces
			if (n0_defined == 1) {
				double ffgamma2_C1 = Ufact[2] * (  ax[an[0]][j]/R - CF7*Re[j] );
				double ffgamma2_C2 = Ufact[2] * (  CF8*ax[an[0]][j] - Re[j]/l[0] );

				ffgamma2[2] =   ffgamma2_C1; 
				ffgamma2[3] =   ffgamma2_C1; 
				ffgamma2[0] = - ffgamma2[2] - ffgamma2_C2; 
				ffgamma2[1] = - ffgamma2[2] + ffgamma2_C2; 
			}

			force_atom[0][j] = ffR + ffr20 + ffr30;
			force_atom[1][j] = ffR + ffr21 + ffr31;
			force_atom[2][j] = -ffR - ffr20 - ffr21;
			force_atom[3][j] = -ffR - ffr30 - ffr31;

			force_atom[an[0]][j]   += ffgamma0[0] + ffgamma1[0] + ffgamma2[0];
			force_atom[1-an[0]][j] += ffgamma0[1] + ffgamma1[1] + ffgamma2[1];
			force_atom[an[1]][j]   += ffgamma0[2] + ffgamma1[2] + ffgamma2[2];
			force_atom[5-an[1]][j] += ffgamma0[3] + ffgamma1[3] + ffgamma2[3];
		}

		for(int atn=0;atn<4;atn++){
			for(int j=0;j<3;j++)
			{
				forces[atn*3+j] = force_atom[atn][j]*0.00138; // forces in N*10^-10
			}
		}

		return PE; // U in K
	}


	void CalcAtomCoordByAngles(double xc[3], double r0, double r1, double R, double g0, double g1, double g2, double *x)
	{

		// 0
        x[0] =xc[0] + r0 * cos(g2);
        x[1] =xc[1] + -r0 * sin(g2);
        x[2] =xc[2];
        // 1
        x[3] =xc[0] + -r0 * cos(g2);
        x[4] =xc[1] + r0 * sin(g2);
        x[5] =xc[2];
        // 2
		x[6] =xc[0] + R - r1 * cos(g1);
        x[7] =xc[1] + r1 * sin(g1) * sin(g0);
        x[8] =xc[2] + r1 * sin(g1) * cos(g0);
        // 3
        x[9] =xc[0] + R + r1 * cos(g1);
		x[10] =xc[1] + -r1 * sin(g1) * sin(g0);
		x[11] =xc[2] + -r1 * sin(g1) * cos(g0);

	}

	void Rotate(double r[3],double w[3],double ang, double *res)
	{
		double cosa =cos(ang), sina = sin(ang);
		res[0] = (cosa+(1-cosa)*w[0]*w[0])*r[0] + ((1-cosa)*w[0]*w[1]-sina*w[2])*r[1] + ((1-cosa)*w[0]*w[2]+sina*w[1])*r[2];
		res[1] = ((1-cosa)*w[0]*w[1]+sina*w[2])*r[0] + (cosa+(1-cosa)*w[1]*w[1])*r[1]+ ((1-cosa)*w[1]*w[2]-sina*w[0])*r[2];
		res[2] = ((1-cosa)*w[0]*w[2]-sina*w[1])*r[0] + ((1-cosa)*w[1]*w[2]+sina*w[0])*r[1] + (cosa+(1-cosa)*w[2]*w[2])*r[2];
	}





};

